Approximation properties of the $q$-sine bases

APPROXIMATION PROPERTIES OF THE g-SINE BASES

LYONELL BOULTON 1 AND GABRIEL LORD 2

Abstract. For q > the eigenfunctions of the non-linear eigenvalue prob-
lem associated to the one-dimensional g-Laplacian are known to form a Riesz
basis of L 2 (0,1). We examine in this paper the approximation properties of
this family of functions and its dual, and establish a non-orthogonal spec-
tral method for the p-Poisson boundary value problem and its corresponding
parabolic time evolution initial value problem with stochastic forcing. The
principal objective of our analysis is the determination of optimal values of q
for which the best approximation is achieved for a given p problem.

Contents

1. Introduction

2. The q-sine basis and its dual

3. Approximation of source terms

4. Numerical solution of the p-Poisson equation

5. The time dependent p-Laplacian
Acknowledgements

References

2

m
m

EH

mi
us
m

Date: February 2011.

1

2

1. Introduction

The p-Laplace operator or p-Laplacian, a generalization of the ordinary Laplace
operator, arises naturally in applications from physics and engineering including:
slow-fast diffusion related to particles [2] , superconductivity [5] , wavelet inpainting
[23] , image processing [17] and game theory [22] ■ A typical application in the large
p limit is a model for slow-fast diffusion for sandpiles [14] . Recently there has been a
significant amount of research activity encompassing methods of approximation for
solution of non- linear partial differential equations involving this operator [HI HI E] -
The aim of the present paper is to further contribute to this activity by considering
the particular case of the one-dimensional p-Laplacian and examine in detail the
approximation properties of a generalized spectral method described as follows.

Let p > 1. For z G K let [z] p_1 = z\z\ p ~ 2 . By extension from the linear
case corresponding to p = 2, we define the one-dimensional p-Laplacian to be
the differential operator A p u = Here u : [0,1] — > K is such that

G i? 1 (0, 1). The corresponding p-Poisson boundary value problem is given

by

A v u(x) = g(x) < x < 1

(1) P ~ ~
y ' u(Q) = u(l) =

where g G £ 2 (0, 1). We also consider the related evolution equation

r A , f „ , dW( X ,t)

— — — = [A p u{x, t) - g(x)\ + v —

( 2 ) u(x,0) = 0<x<l
u(0,t) = u(l,t) = t>0

that includes a stochastic forcing term where the noise intensity v > and W is a
space-time Wiener process. Here

(3) W(x,t) = J2l 3 n(t)Mx)

nGN

where /3 n are independent scalar Brownian motions, ip n are the eigenfunctions of
the covariance operator Q and a n the corresponding eigenvalues, see [101 12T1 [T^] .
For the case of space-time white noise we have Q — I. In the case of a deterministic
system when v = 0, ([I]) is the steady state solution of

Let q > 1. The so called q-sine functions are defined as the eigenfunctions of
the g-Laplacian eigenvalue equation [131 023 HO] :

-Mx) = ( g -1)A[ U ( 3 :)]'- 1 0<x<l
u(0) = u(l) = 0.

This family of functions and the corresponding problem Q was studied over 30 years
ago by Elbert P~3] and later by Otani [5U], Bennewitz and Saito [7] in connexion
with the computation of optimal constants in Sobolev-type inequalities. They gen-
eralize in a natural fashion the 2-sine basis (corresponding to the linear case), and
they have very similar periodicity and interlacing structures.

In [8] analogues of the classical completeness and expansion theorems for the
g-sine functions were established for q > 12/11. Specifically it was shown that they
form a Riesz basis of L 2 (0, 1). This leads to the following question: what are the

3

approximation properties of this basis, as well as its dual basis, and how they relate
to the approximation properties of the standard 2-sine basis?

Below we address this question by examining approximation of the solutions of
(JlJ and ([2| via projection methods with a g-sine and a dual g-sine basis, regarding
q > 1 and p > 1 as free parameters. A main focus of attention is the determination
of optimal values of q for which the highest order of convergence is achieved in
a given p-problem. We demonstrate that standard properties of the 2-sine basis
applied to the p — 2 problem (such as super-polynomial convergence when g(x) is
smooth) are lost, when a g-sine basis for q ^ 2 is considered for ap^2 problem.
As it turns out, the property of being a basis for the g-sine functions conceals a
remarkably rich structure which is far from evident given the apparent simplicity
of problem Q.

Background material on the p-sme basis and its dual is considered in ££2j There
we examine a matrix representation of the Schauder transform introduced in [5].
This will be crucial in our subsequent analysis as it gives rise to a stable procedure
for constructing numerically both bases.

In <|3] we find estimates for the approximation of square integrable functions in
terms of their regularity. The dual g-sine basis turns out to have very similar ap-
proximation properties as the q = 2 basis (lemma |2]). On the other hand, however,
it is fairly simple to construct smooth functions such that their g-sine Fourier coef-
ficients do not decay faster than a power —5/2 for q > 2 (lemma [3]). The latter is
in stark contrast with the most elementary results in the numerical approximation
of solutions of differential equations by orthogonal spectral methods.

Section [4] is devoted to the p-Poisson boundary value problem. In theorem [5] we
find explicit uniform bounds on the distance between any two solutions of ([I]), given
the distance between the corresponding right hand sides. We then examine in detail
the numerical computation of solutions of for source terms that are subject to
various different regularity constraints. As it turns out, the estimates established in
Theorem[5] appear to be sub-optimal. A more thorough investigation in this respect
will be reported elsewhere. See [TT] for related results in the context of finite element
approximation of the solutions of ([!]), including the higher dimensional case.

In the final ^5] we study the numerical approximation of solutions to ([2| both in
the deterministic and stochastic systems. We describe our discretization strategy
and solve this problem for different values of p and v. Our results provide evidence
on the performance of the g-sine basis for the solution of |2]), by showing the
dependence on the parameter q of numerically computed L 2 residuals.

2. The g-siNE basis and its dual

The g-Laplacian eigenvalue problem Q, although non-linear, has a fairly simple
structure. The eigenvalues are found to be A = {nir q ) q where ir q = qsi ^/ q ) ■ The
first eigenfunction fi(x) associated to the first eigenvalue (ir q ) q is strictly increasing
in [0, 1/2], decreasing in [1/2, 1] and it is even with respect to x = 1/2. It can be
extended to an odd function (with respect to x = 0) in the interval [—1,1] and then
to a 2-periodic C 1 function of R. If q > 2 then f"(x) is singular at x = 1/2. The
eigenfunctions f n (x) associated to the eigenvalues (nir q ) q satisfy f n (x) = fi(nx)
for all n > 2.

In figure [I] we have plotted the first three eigenfunctions (top) and their corre-
sponding derivatives (bottom) for (a) q = 1.4 and (b) q = 10. They typify the

4

case q < 2 (a) and q > 2 (b), respectively. For large q the basis functions /„ ap-
proach zig-zag functions, which are the eigenfunctions of the oo-Laplace eigenvalue
problem.

Below we always assume that the family {f n }neN is normalized by the condition
/4(0) = nn q and leave implicit the dependence of f n on q. In the special case q = 2
we write e n (x) — \/2 sin(n7ra), so that {e„}„ e N is an orthonormal basis.

(a)

(b)

— J,

— V

V: \

(0

t •

>

1 -5

S

■o

-10

Figure 1. Approximation of the q-sine functions fj for j = 1,2,3
(top) along with their derivatives (bottom) for (a) q = 1.4 and (b)
q = 10.

The Pythagorean identity generalizes to the q-sine functions [T3] as
(5) [/i(ar)| fl + 7r-«]/i(a:)|« = 1.

Integrating this differential expression for small enough x leads to the following
explicit representation for the inverse function

As we will see below, this representation plays a crucial role in the numerical esti-
mation of f n (x).

Let the Schauder transform, T q , be the linear extension of the mapping e„ i — >
f n . Then T q : £ 2 (0, 1) — > L 2 (0, 1) is an invertible bounded operator for all q >
12/11, [U Theorem 1]. Thus {f n }neN is a Riesz basis of L 2 (Q, 1) for such range of
the parameter q. Further evidence presented in [5] suggests that in fact this is also
the case for all q > 1, but at present this has not been proved rigorously. Unless
otherwise specified we will assume from now on that q > 12/11.

The property of a Riesz basis ensures that every g e L 2 (0, 1) is represented by
a unique series expansion g = X^^Li a nfn which is convergent in norm. The q-sine
Fourier coefficients, a n £ K, are given explicitly by a n = (g, /*) where {/^}neN is

, 3 ,T*f*) for all

the basis dual to {/„}„ g n- Since S jk = (fj,f%) = (T q ej,f%) = {
j, k e N, then /* = (T~ 1 )*e n . It turns out that /* 7^ /„ for q ^ 2.

The following matrix representation of T q is fundamental to our analysis. Let

TqU) = AO')

v2 / f\ (x) sin(j7ra;) da;

5

be the jth 2-sine Fourier coefficient of f\{x). Then the fcth 2-sine Fourier coefficient
of f n (x) is given by

fn(k) = V% / fi(nx) sin(fc7ra;) da;
Jo

(7) =2 fiifn) / sin(m7rna;) sin(fc7ra;) dx

m odd ^°

{T q (m) if mrt = fc for some m odd
otherwise.

Hence T q e n = Xim=i T ?( m ) e mn an d therefore T q has a lower triangular matrix
representation in the orthonormal basis {e n } ne jq. See figure [2]- (a).

Figure 2. In (a) we plot the distribution of the non-zero entries
of a 1000 x 1000 truncation of T q . The insert corresponds to a
20 x 20 truncation. The matrix entries are constant along each of
the "quasi-diagonals" seen in the picture. In (b) we plot s against
q where s is the numerically estimated H s regularity of the basis
function fi{x). Note that at q = 2 the regularity is infinite.

Remark 1. The basis of eigenvectors of the oo-Laplace eigenvalue problem are
zig-zag functions, [HI Section 5]. In this case we can write ^(j) explicitly. As it
turns out,

g

lim r q (j) = = Teoij).

q— >oo 7 7T

Let s > 0. Below we denote by iJp er [0, 1] the Sobolev space of 1-periodic func-
tions g £ i 2 (0,l) such that £^(1 + j 2 ) s \g(j)\ 2 < oo.

Lemma 1. Let f\ be the first q-sine function as defined above. If 1 < q < 2, then
fl G H p 2 cr [0,l]. Ifq>2, then f t £ (J s <3/2 #per[0, 1] \ ff p 3 cr [0, 1] • Ifq> 4, then
h $ H* el [0,l}.

Proof. A straightforward argument involving integration by parts yields

2y/~2 f 1 / 2

Tqii) = — f"(x)sin(j<7rx)dx.

3 k

6

Then

(8)

3 71 Jo 3 77 Jo

_ 2y2 , 1/2 _ 2V2ir q

„'2_2 Un x )\o — -2„2 •

1/2 |/r(x)| 2 dx<l max \h(y)f = 4:

Hence r q {j) = o(j~ 2 ) as j — > oo and /i 6 U s<3 ^ 2 iJp Cr [0, 1]. From |5| it follows

that /('(z) - ft(/i(a;)) for h(y) = -v^V^l - y")^ .
Let 1 < g < 2. Then

1/2

lf,"WI 2 rlT< ......

2 ye [0,1]

so A ei^ er [o,i].

Let q > 2. Since lim^ii /"(x) = Too, then f x g i?p Or [0, 1]. Moreover, f x {x) >
2x for < x < §. Then

l/"W| 2 dx > j 1 ' 2 {2x) 2 ^ 1 \l - (2a;) 9 ) ^ dx.

If g > 4, the integral on the right hand side diverges and so /i ^ -ffp CI .[0, 1]. □

Evidently the H s regularity for f\ (x) found in lemma [I] is not optimal. Fig-
ure |2j-(b) shows a numerical estimation of the precise value of s(g), such that
fi £ iJp Or [0,l] for r < s(q) and fi g Hp CI [0, 1] for r > s(q). The data for this
graph was obtained by computing the decay rate of T q (j) for a large truncation of
T q . A thorough investigation closely related to this lemma in the higher dimensional
context can be found in [T^] and references therein.

In the large q we have a limit of s = 3/2 and this is confirmed by Remark [I] At
q = 2 we simply have f± = sin(7rx) and so f\ is in ifp Cr [0, 1] for all s. According
to lemma [l] the curve should remain below s = 3 as q — > 2 + . For q > 3 the graph
suggests /i ^ Hp C1 .[0,l]. However, if q < 3 it suggests s(q) > 2. As q — > 1 the
regularity drops and the limit seems to approach s = 2. There is an interesting
"peak" of regularity around q = 1.6 which we can not presently explain.

The 2-sine Fourier coefficients of /^(x) are given by the n-th columns of (T^ 1 )*.
This operator has an upper triangular representation in the 2-sine basis. Then
/*(x) are trigonometric polynomials of order n. In fact, /|(x) = r ? (l) _1 sin(27rx)
are parallel for all q > 1. Unlike the g-sine functions, not all dual g-sine functions
have the same periodicity structure.

We now describe a stable numerical procedure for computing these two bases.
The g-sine functions can be approximated by first estimating /i using Numer-
ical integration yields /-f 1 on [0, 1/2]. Although the integral is singular at y = 1,
the value /^~ 1 (1) = 1/2 is known, so we do not need to consider quadrature points
too close to this singularity. In our numerical procedure we chose a fine uniform
grid and apply a cumulative Simpson's rule. This gives an approximation fx of
/i for x £ [0,0.5] with a controlled tolerance and by symmetry we obtain f\ for
x € [0, 1]. Note that /i is given on a non-uniform grid. The Pythagorean identity
([5]) immediate yields the derivative /{ of f\ and hence we can use it to approximate
the former with a /{, defined also on the non-uniform grid.

7

Once we have constructed /i and f[, we use periodicity and symmetry to find
corresponding approximations of /„ and f n for n > 1 . This involves considering
scaled copies /i and /{, to form /„ and j' n . Here the number of non-uniform grid
points on [0, 1] grows with each n = 2 : N.

To obtain /„ and f n on a uniform grid Xj — jh for j = : J, rather than the
non-uniform grid that arises from the numerical integration, we have considered
numerical interpolation by piecewise cubic polynomials. This gives {fn}n=i an( f
{(fn) h }n=i defined at x — Xj. The former is the approximated basis and the latter
the corresponding derivatives that we use for further computation.

Figure [I] was generated with an implementation of the numerical scheme just
described on an uniform grid with J = 4000 points (h — 2.5 x 10~ 4 ). The integral
(|6| was approximated with 2 x 10 5 points.

The dual basis is found from an N x N truncation, , of the Schauder trans-
form. In practice, we first compute T q (j) and assemble . Then we define approx-
imations (fn) h for n = 1 : iV, as the trigonometric polynomials whose kth 2-sine
Fourier coefficients are the (n, k) entry of the matrix (T^) -1 .

Remark 2. In our numerical approximation of the set of basis functions we are
careful to fully resolve oscillations on the basis function fx, taking at least 20 mesh
points per wavelength and for most computations 100 per wavelength. For N < 50
we use h = (1007V)" 1 and so resolve each oscillation in the basis function with 100
spatial points. For N > 50 we use h = (20N)~ 1 and so resolve with 20 points.
We also examined convergence of orthogonality of the basis and dual in the spatial
discretization and noted 0{h 2 ) for q « 10 through to 0(h 4 ) for q 1.4.

3. Approximation of source terms
Any given g £ L 2 (0, 1) can be approximated by either

JV N

9(x) « g* N (x) = ^2(g,fj)fj(x) or g(x) « g N (x) = ^(g, /*) fj(x)

for large N. These two expansions converge as N — > oo in the norm of £ 2 (0, 1) and
also pointwise for almost all x 6 [0, 1]. Unlike in the linear case corresponding to
q = 2, the rates of decrease of || g — g^\\ and \\g — gj^W can be very different when
q^2.

Since the dual basis {/^} n eN comprises trigonometric polynomials, on the one
hand we can formulate the following natural statement.

Lemma 2. Let g e L 2 (Q, 1). For all JVeN,

\9-9n\

<

2V2

N

9(n)e n

Proof. By definition g — g^ — X)^Ljv+i(^?*5' e «)(^ 1 g) le «- Since the matrix asso-
ciated to T* is upper triangular (see Section [2l,

(T*g,e n ) = J2[T*] nk g(k) = T q (k)g(nk) .

k=l

8

According to Q, we have Y^k=i l T ?(^)l — ^73- Hence

l(T;5,e„)| 2 = |5> g (fc)5(nfc)| = I £ t,(*0 1/2 t 9 (*) 1/2 ?(**)|

k=l k=l

oo

<

(E M fc )i) (E W*)iis(»*)i a ) < ^= E w*)ii£(»*)i s

fc=l fe=l V fc=l

Thus,

OO OO

n=JV+l v n =iV+l fe=l

00 00

K'w'vkY. E i^wii^i^ii^ii^Em^i E i^)i s

2 fc=ln=JV+l Zy/Z k=l n=N+l

2 00

2

<II^TA£k(*)l E I^HI^Tf E Wi\

V k=l n=JV+l 71=^+1

□

Therefore the q-sine dual expansion of any g <= C°° converges super-polynomially
fast. On the other hand, however, it is not difficult to construct examples of smooth
functions g with a subsequence of q-sine Fourier coefficients decaying slowly.

Lemma 3. Let g{x) = e±(x). If j is prime, then (g,f*) — T g (j)r g (l) .

Proof. Since T q has an (infinite) lower triangular matrix representation in the or-
thonormal basis {e n } n6 N, we can find the entries of the corresponding matrix rep-
resentation of T^ 1 by pivoting and forward substitution (Gaussian elimination), ft
is readily seen that T^ 1 is necessarily lower triangular and its diagonal should be
constant and equal to r q (l)~ l .

Assume that j is prime. According to ([7]), the only non-zero entries in the jth
row of T q are r q (j) in the first position and r g (l) in the jth position. Therefore,
the only non-zero entries in the jth row of (T g ) _1 are — in the first position
and r 9 (f) _1 in the jth position. As the Fourier sine coefficients of f* are obtained
from the jth column of (T*) _1 , the desired conclusion follows. □

By virtue of lemma [I] the prime g-sine Fourier coefficients of sin(7rx) for q > 2
can not decrease faster than j~ 5 ^ 2 in the large j limit. Observe that this is in
stark contrast with the most elementary results in the numerical approximation of
solutions of differential equations by orthogonal spectral methods.

Remark 3. The finite set of basis functions /„ and dual /* for n = 1 : N generate
corresponding N dimensional subspaces Vat and of L 2 (0, 1). Instead of com-
puting directly with these non-orthogonal bases, one can apply the Gram-Schmidt
algorithm in order to obtain orthonormal bases of these subspaces. This has a
numerical advantage of not needing to store both the basis and dual. We also
considered this approach, however we found little advantage in terms of accuracy.

Let us now consider various numerical tests on the approximation of regular
functions by {/ n }neN and {/^} n eN- Once the bases have been obtained on uniform

9

mesh, we can examine their approximation properties numerically by looking at
the decay of suitable residual for benchmark sources g € L 2 (0, 1).

Figure [3] shows the typical outcomes of an experiment to determine the depen-
dence in q of the L 2 (0, 1) residual. We have fixed here N = 40, varied q = 1 : 100,
and computed residuals in the approximation of the following four functions:

Sa(a0 = /i(a0 + (2.5)/io(a0 for q = 10,
9b(x)

x e [1/4,3/4]
otherwise

(9)

=0*0

(7/3)x x € [0,3/7]

-{2l/4)x + 13/4 x e [3/7, 7/14]

(21/4).t-2 x € [7/14,4/7]

-a; + 11/7 x €[4/7,1]

g d (x) = (q-l)(kn q rf 1 \f 1 \^ for q = 3,

in (a)-(d) respectively. In the figure we include both the basis and its dual, as well
analogous calculations with the orthogonalized bases of V40 and V^ . In (a) we see
an optimal q op t = 10 as expected, in (b), (c) and (d) we increase the regularity of g
and see g opt decrease with values of 4.25, 2.9, 2.55 for (b), (c) and (d) respectively.

This experiment gives a general insight about the g-behavior of L 2 -residuals in
the approximation of functions with different degrees of regularity by g-sine bases
and their duals. For the simple functions considered, our calculations indicate the
following general behavior of the g-sine basis:

(i) as q — >• 1 the residual deteriorates,

(ii) there is always a single minimum corresponding to an optimal q = (fopt,

(iii) as q —¥ 00 the residual curve becomes asymptotically constant with no local
maximum for q > 2.

In contrast, the approximation error is almost constant in q for the g-sine dual basis.
This is indeed a consequence of Lemma [2] and the fact that /* are trigonometric
polynomials. Our tests indicate that there does not seem to be a clear advantage
in orthogonalizing the basis or the dual basis for errors with an order of magnitude
above 10 -5 .

g{x)

9b{x)
9c{x)
9d(x)
sin(7rcc)
Table 1.

g = 1.8
-0.4862
-1.4266
-1.9368
-2.2778

<7 = 2 (exact)
-0.4905 (-0.5)
-1.4474 (-1.5)
-1.9952 (-2)
NaN (-00)

9 = 3

-0.5005

-1.4959

-1.9826

-1.9560

9 = 5

-0.5052

-1.4842

-1.5934

-1.6237

g = 10

-0.5069

-1.4438

-1.4595

-1.4988

Rates of convergence in N for g with different degrees

of regularity and selected values of q.

In table[T]we have estimated a > such that ||<? — g^\\ < (3N~ a where /3 > is
independent of N. The data indicates that for g £ Hp CI [0, 1], q pt ~ 2. Moreover,
as q increases, we should expect a to decrease and stabilize always below 1.5 for
g(x) = sin(7ra;). This is indeed suggested by lemma [I] figure [2]-(b) and lemma [3]
and it is confirmed by the last row of the table.

In figure[4]-(a) we have computed the residual \\g — <7io|| (red) and \\g — 520II (blue)
for randomly generated g € ifp C1 .[0, 1] constrained to g(0) — g(l) = 0. To construct

10

(a) (b)

Figure 3. We vary a g-sine basis and dual basis and examine how
the residual changes using N — 40 modes for (a) g a combination
of two 10-sine basis elements, (b) a piece-wise constant function
<Zopt ~ 4.25, (c) a piece-wise linear continuous g pt = 2.9 and (d)
a differentiable with discontinuous derivative function q opt — 2.55.
See ^.

a random function g such that the i/p OI .[0, 1] norm is finite, we find g~(j) — ajfij
where a,j = (1 + j 2 )~ s ^ 2 \j\~^~ s for some small 5 > and f3j are independent
identically distributed N(0, 1). Two realizations of functions g obtained in this way
are shown in the inserts in figure [4j Note that q opt is achieved at different places
but close to 2. In figure |l]-(b) we examine q op t over 200 realizations of functions in
ifp O1 .[0, 1]. For any fixed realization and value of N, q opt ^ 2 although in the limit
this optimal parameter appears to be close to 2. For N = 100 we also show the
results of 1000 realizations. In this case, the mean value is closer to 2 although the
variation is still large. In figure |4]-(c) we show the mean values where functions are
taken in H poT [0, 1] for s = —0.5, 0, 0.5, 1 and 2 with 200 realizations. For s < 2 the
variability in q opt is far larger. For fixed N, q op t > 2.

The observed outcome of this experiment strongly support the conjecture that
<fo P t typically approaches 2 as the regularity of g is increased.

11

(a) (b) (c)

Figure 4. In (a) 1? residual as a function of q where a g-sine
basis is used to approximate two different randomly generated g £
Hp CI [0, 1] satisfying g(0) — g(l) = 0. Note that the optimal g op t
occurs at different values. In (b) and (c) we examine q opt as N
increases. We include the mean and standard deviation over 200
realizations of random functions subject to the same constraint.
For N — 100 we also show the results from 1000 realizations.

4. Numerical solution of the p-Poisson equation

We now address the question of approximating the solutions of ([I]) by means
of a g-sine and dual basis. In view of lemmas [2] and [3j we begin by determining
uniform estimates on how sensitive this solution is under perturbations of the right
hand side. Analogous questions have certainly been considered in more general
frameworks, however here we focus on the explicit calculation of the constants
involved.

A key ingredient in the estimates presented below is the fact that can be
integrated explicitly. Let the Volterra operator

Vg(x) = [ g(t) dt.
Jo

Note that Vg £ H 1 ^, 1) for all g £ L x (0, 1). Furthermore V : L s {0, 1) — > L*(0, 1)
is a contraction operator for all 1 < s, t < oo and its norm can be explicitly
determined [71 Theorem 1.1]. Let

Jo

Then hgfa) is a continuous function, decreasing in 7, for all fixed g £ L (0, 1). Let
min Vg(x) < 70(g) < max Vg(x)

x£[0,l] x£[0,l]

be the unique root such that /i s (7o(s0) = 0. Then

(10) u(x) = J lVg(r) - 70(5)]^ dr = V (lVg{r) - 7 o(<?)]^) (*)

is the unique solution of ([!]).

We firstly establish concrete Holder estimates on h g {^). Without further mention
in this section we will fix r = and denote by || • || s the norm of L s (0, 1). In the
case s = 2 we will continue suppressing the sub-index.

12

Lemma 4. Let g <E L 1 (Q, 1) and — \\g\\i < 7 < [i < ||.9||i- TTien

(M-7)< l f^[h g (j) - h g (^)} 0<r<l

(M-7) r <2'- 1 [^(7)-^(/i)] r>l.

Proof. Suppose first that < r < 1. From the graph of Jz] r for \z\ < M it is
readily seen that [z] r - (wf > rM r - x (z - w) for all -M <w < z < M. Then

(/i - 7) = / [(V 5 (r) - 7) - (Vs(r) - /x)] dr
Jo

^ — — M7) - KM]

for M = 2||g||i, and 7 and /1 as in the hypothesis.

In a similar fashion, let r > 1. Then (z — w) r < 2 r ~ 1 ([[,z] r — Jw] r ) for all
—M < w < z < M. Indeed, if < w < z, a straightforward argument shows that

> z T ~ x > (z-wy- 1 ;

.6 — OJ ^ — UJ

if w < z < 0,

and if w < < z,

> M > (z- wY

z — w

■ HI - MT • z r ~\w\ r 1

mm — ; r — = mm ■ —

T>6 (z~w) r T>6(z+\w\) r 2 r " 1
achieved when w = — z. Thus

(/x - 7 )' r = /' [(Vg(r) - 7) - (Vff(r) - M )f dr

< T- 1 [^(7) - ■

□

Theorem 5. Let u and u be solutions of ([T]) with corresponding sources g and g.
Let m = max{||g||i, ||<?||i}. Then

\\u - u\\ < 2'- r (\\g - g\\ + (4 "^HIs -g\\l) < r < 1

||u - fi|| < r2 r - 1 m r - 1 (\\g - g\\ + 2 2 -^ r m 1 - 1 / r r 1 / r \\g - g\\{ /r ^

r > 1.

Proof. Let < r < 1 and s = 2/r > 2. By virtue of ([10 1

||«-fi|| < ||^([^-7o(g)] r -[^"7o(ff)F)|| s < \\lVg-jo(gW-iv~g-io(~gW\

Note that \{zj r - {w} r \ < 2 1 ~ r \z - w\ r for all < \w\ < \z\. Thus

Hl^ff " lo(gW ~ IV g - 7o(ff)]1 s 1/r < 2^ (\\V(g - g)\\ + \jo(g) - 7o(<?)|)

<2^ r (||.g-.g|| + |7 (5)- 7o ( ff )|).

13

According to lemma [2J

ho(~g) - 7o(ff)l < l^(7o(5)) - h g ( l0 (g))\

(2m) 1 "

(2m)

<

(4m) 1 -'' „

< — — lis-

[^-7o(5)] r -I^-7o
(4m) 1 ~ r

1^(70(9)) -h- g (^(g))\

dr

g\\ r r <

~g\\ r v

r r
This ensures the first statement.

Let r > 1. In a similar fashion as before, we see that

||«-fi|| < ||I^ 5 -7o(.9)] r -[^-7o(5)]1

< r2 r ' 1 m r - 1 (\\V(g - g)\\ + \ l0 (g) - lo (g)\)

<r2 r - 1 m r - 1 (\\g-g\\ + \ lQ (g)-lo(~g)\)-
Lemma [4] and similar arguments as for the previous case, yield

bto® - 7o(ff)l r < r- 1 \h- g { l0 {g)) - h g ( l0 (g))\

<2 2r - 2 m r - 1 r\\ g - g \\ 1 .

This completes the proof.

□

The right hand side bound in the above theorem approaches 2 as r — > indepen-
dently of the value of \\g — g\\. In fact \\u — u\\ < ||Mli f° r K — Wg~ 7o(ff)] r ~ W(3~
7o(g)] r . Since l r (r) -> for almost all r £ [0,1] and |Z r (r)| < 2 + \\Vg\\ 00 + \\Vg\\ 0C ,
the Dominated Convergence Theorem yields \\u — {t|| — > as r — > 0. This is a well-
known property of the oo-Laplacian, see |22j . As mentioned in the introduction,
the approach considered in in the context of finite element approximation of
the solutions of ([!]), may provide an insight on whether the constants found in the
above theorem are optimal.

(a)

10

q - basis

2 10

(b)

Figure 5. Solving the p-Laplacian problem with p = 5. (a) The
most accurate basis for a solution u(x) — sin(7ra;) is the standard
2-sine basis, (b) However for a solution u(x) = fi(x) with q = 5
the 5-sine basis is the most accurate.

We now describe how the p-Poisson problem is discretized. The strong formula-
tion ([T]) leads to the following weak formulation using integration by parts and the

14

boundary conditions:

(11) f |w'| p "Vv'd:r = f gvdx

Jo Jo

for any absolutely continuous test function v. We expand both u and v using a basis
{4> n }n£Th where <p n is either /„ or /*. After truncation this leads to the following
nonlinear system of equations for unknown coefficients Cj = (u, <j)*)\

p-2

M'j^ = (9>^) j = l:N.

The case p = 2 evidently reduces to the standard linear system to solve for the
Poisson equation.

(a) (b)

10 1 10 Z 10 1

q - basis q - basis

Figure 6. We now examine a piecewise constant / for different p
(N = 40). In (a) p = 1.8 and g Q pt = 2 (b) p = 3 and (fopt ~ 2 with
another minimum at q = 2.95, in (c) p = 5 and q opt ~ 3.9 and in
(d) p — 10 and q op t ~ 5.75.

Numerically the right hand side of (12) is approximated via quadrature rules.
The nonlinear system may then be solved for example by Newton's method. Since
the Jacobian is a full matrix in this banded approximation appears to

provide sufficient accuracy. However, the results presented below were found using

15

(a) (b)

10 1 10 1
q - basis q - basis

Figure 7. Comparison of N = 20, 40, 60, 80 in (a) for p = 5 and
g = 1 (b) for p = 10 and g = 1 (also including N — 100). In (c)
we show for g = g^ and p = 5, and in (d) for g = g^ and p = 10.
We note that for N sufficiently large q opt = 2, and a more accurate
solution may be found for smaller N at q op t ^ 2.

the trust-region dogleg method with a full Jacobian matrix as implemented in
Matlab.

In figure [5] we have use this scheme to solve the p-Laplacian problem with p — 5
and examine the L 2 error taking N = 40. In figure [El (a) we have fixed a solution
it(:c) = sin(7rx) where we clearly expect and observe that q = 2 is the optimal basis.
In figure [51(b) we have fixed u(x) = fi(x) for p — 5 and so we observe the q — 5 as
the optimal basis.

In figure[6]we take g — gb from ([£]) which turns out to be a typical form of forcing
for sandpile problems. We solve the p-Laplacian problem for (a) p — 1.8, (b) p = 3,
(c) p = 5 and p = 10 with N = 40 and examine the L 2 error in the solution as we
vary the q basis. To estimate this error we take as exact the solution with 2-sines
and 2N modes. We observe that the optimal basis for representing the solutions is
no longer the standard q — 2 for p — 3, 5, 10. For problems with moderate p (for
example p = 3) we see two distinct minima. For larger p problems however, the
q = 2 basis becomes less competitive.

In table [2] we give estimates of q op t for g — gh and g — 1 for N — 20, 40, 60 and
N = 80 modes. The changes in <7 op t with N are explained by the interchange of the
two minima in figure [6l For N — 100 this interchange occurs for p = 5 and p = 10
as illustrated in figure|7l(b).

16

g

p= 1.8

p = 'S

p = 5

p = 10

<?b (N = 20)

2.05

2.95

3.9

8.0

g h (N = 40)

2.05

2.0 (2.95)

3.9

5.75

ffb (N = 60)

2.2

3.0

4.0

5.95

<?b (N = 80)

2.05

2.0

3.9

5.60

g = l(N = 20)

2.0

2.3

3.1

6.43

.g = 1 (AT = 40)

2.0

2.0

3.0

4.85

g = l(N = 60)

2.0

2.0

3.0

4.7

g = l(N = 80)

2.0

2.0

2.0 (2.95)

4.7

Table 2. Optimal q opt for four different p-Laplacian problems
with either a g — gy, from (|9j) or g = 1. We give estimates for
different values of N. The change in q opt as iV increases occurs as
the two minima interchange, see figures [6] and [7j

5. The time dependent p-Laplacian

In this final section we consider the evolution equation ^ both in the deter-
ministic [y = 0) and the stochastically forced regime (i/ ^ 0). Slow- fast diffusion
is often taken in the large p limit, as a model for sandpile growth. For further
details including arguments about the validity of the modeling see for example
[TBI [H [TBI HH HJ [5] • Our purpose here is to consider the time dependent problem
for different values of p and examine the choice of basis in the formation of the
sandpile.

We discretize the weak form of the equation and truncate to solve the time-

dependent version of ( 12 1 given by

(13)

dt

N

E

fc=i

Cfc

N

E

1=1

C((j

p-2

(b'k^dx - (g, fy) + v(^A 3 ) j = 1 : N.

This expression is then discretized in time to get a nonlinear system of equations
to solve for c™ at each step. Below we consider an Euler's method which reduces to

„n+l

At

N

E'

fc=i

N

E

i=i

,n+l.

P-2

M

E

bj dx

1/2

(x)A/3„

where A/3 m are independent identically distributed random variables with mean
zero and variance At (recall ([3])), and we have assumed the eigenfunctions of Q are
now given by e n . For the case of stochastic forcing we take M » N.

In figure [8] we consider v = with g = g^ and we plot the time evolution for
(a) p - 1.8 (<fo p t = 2), (b) p = 3 (cfc p t = 2.95), (c) p = 5 (<z opt = 3.9) and
(d) p = 10 (<7opt = 5.75). In each case the evolution found numerically quickly
converges towards the steady state.

Figure [8] was obtained by choosing in each case the predicted q = q op t basis
for the steady state. Intuitively this choice of basis should be near optimal for
large t. On the other hand however, there is no reason to presume that it is so

17

(a) (b)

Figure 8. Solution of ^ with g = g^ for (a) p = 1.8 using q = 2,
(b) p = 3 using q = 2.95, (c) p = 5 using g = 3.9 and (d) p = 10
using g = 5.75. As time evolves the solutions approach the steady
state which, for increasing p, turns out to be close to a hat function.

(a) (b)

Figure 9. Residual in time for the numerical solution of We
choose v = 0, N = 40 and take as a true solution one computed
with N — 100 and q = 2. We fix p = 10 and compare basis with
q = 1.8, q = 2, q = q opt and q = 10. Here (a) corresponds to
9 = 9b (q op t = 5.75) and (b) to g = 1 (q opt = 4.85). We see
that asymptotically in time g opt has the smallest error. At a few
points in time during the transient state, the q = 1.8 basis is more
accurate than the q = 2.

for small values of t. In figure [9] we examine the spatial error at each step for
four different bases with N — 40. For comparison we have chosen p = 10 and
q G {1-8, 2, q opt , 10}. As the solution approaches the steady state, in each case the
basis q pt has the smallest error as is expected (in fact a gain of almost an order of
magnitude in accuracy compared to q = 2 is observed). It is remarkable however
that, in the transient regime, these bases appear to be far less accurate than other
choices of q, such as q = 2 and q — 1.8.

Let us now consider the stochastically forced case, v > 0. Figure[l0]shows plotted
solution with time dependent noises which are: H 1 in space and white in time for

18

(d)

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

(e)

(f)

2 0.4 0.6 0.8 1

Figure 10. Sample realization with p = 10 of the time evolution
problem with a time-dependent stochastic forcing v = 0.2. The
time evolution is shown with equal steps of At = 0.01. In (a)-(c)
we show the effect of an H 1 noise in space, white in time. In (d)-(f)
we show the effect of a white noise in both space and time.

(a)-(c), and white both in space and time for (d)-(f). The forcing corresponds
to v = 0.2. We expect that on average the noisy solution is simply that of the
deterministic system (which has a unique stable solution). We see in (a)-(b) and
(d)-(e) the time evolution of the solution for one particular realization of the noise.
This should be compared with the deterministic case in figure [8j (d).

In figure 10 -(c) and figure 10 -(f) we plot the evolution of the error for fixed time-
step with q £ {1, 2, q op t, 10}. Similar to the deterministic case, in the transient
regime the error for <7 opt and q = 10 is far higher than that for q — 2 or q = 1.8.
Where the noise effects are large (such as for the spatially correlated noise in H 1

in figure 10 (c)), the optimal basis is no longer clear. However where the effect is
small (such as the white noise in figure[lO[(f)), we clearly see in the time dependent
evolution about the steady state that the q — 4.85 outperforms the other bases (and
in particular q = 2). For sandpile-type problems the introduction of a spatially
white noise is a natural choice. It is interesting to note the small effect on the
dynamics in this realization.

At present it is unclear whether the use of a g-sine basis with q =/= 2 provides any
real computational advantage over the natural choice q = 2 for the solution of
There is clearly a computational overhead in obtaining the former for q ^ 2 that
impacts significantly on efficiency. However we have observed that the nonlinear
problem (12 1 is solved faster in the optimal basis. This is certainly worth further
investigation. For example for TV = 40 we observe approximately a 20% speed up
in the nonlinear solve over the standard 2-sine basis. Thus, if solving many fixed
p-problems, the corresponding optimal g-sine basis may not only be more accurate
but also more efficient. In our current implementation the g-sinc basis may be

19

precomputed and stored. Lemma [T] give some prelimeinary indication on how to
solve the problem of apriori determining the optimal basis. Our numerical results
suggest that for large p we expect an optimal basis with q > 2 for a right hand side
with a discontinuous derivative.

Acknowledgements

We kindly thank Adrien Vignes and Bryan Rynne for their thoughtful comments
and involvement in discussions related to this paper. The first author acknowledges
support from the Universite Paris Dauphine, where part of this research was carried
out.

References

[1] F. Andreu, J. M. Mazon, J. D. Rossi, and J. Toledo. The limit as p — > oo in a
nonlocal p-Laplacian evolution equation: a nonlocal approximation of a model
for sandpiles. Calc. Var. Partial Differential Equations, 35(3):279-316, 2009.

[2] G. Aronsson, L. C. Evans, and Y. Wu. Fast /slow diffusion and growing sand-
piles. J. Differential Equations, 131(2):304-335, 1996.

[3] J. W. Barrett and W. B. Liu. Finite element approximation of the p-Laplacian.
Math. Comp., 61(204):523-537, 1993.

[4] J. W. Barrett and W. B. Liu. Finite element approximation of the parabolic
p-Laplacian. SI AM J. Numer. Anal, 31(2):413-428, 1994.

[5] J. W. Barrett and W. B. Liu. Higher order regularity for the solution of some
nonlinear degenerate elliptic equations. SIAM J. Math. Anal., 24:1522 - 1536,
1993.

[6] J. W. Barrett and L. Prigozhin. Bean's critical-state model as the p — > oo limit
of an evolutionary p-Laplacian equation. Nonlinear Analysis, 42(6):977-993,
2000.

[7] C. Bennewitz and Y. Saito. Approximation numbers of Sobolev embedding
operators on an interval. J. London Math. Soc, 70:244-260, 2004.

[8] P. Binding, L. Boulton, J. Cepicka, P. Drabek, and P. Girg. Basis properties of
eigenfunctions of the p-Laplacian. Proc. Amer. Math. Soc., 134(12):3487-3494
(electronic), 2006.

[9] A. Caboussat and R. Glowinski. A numerical method for a non-smooth

advection-diffusion problem arising in sand mechanics. Commun. Pure Appl.

Anal, 8(1):161-178, 2009.
[10] G. Da Prato and J. Zabczyk. Stochastic Equations in Infinite Dimensions,

volume 44 of Encyclopedia of Mathematics and its Applications. Cambridge

University Press, Cambridge, 1992.
[11] C. Ebmeyer and W. B. Liu. Quasi-norm interpolation error estimates for the

piecewise linear finite element approximation of p-Laplace equations. Numer.

Math., 100:233-258, 2005.
[12] C. Ebmeyer, W. B. Liu, and M. Steinhauer. Global regularity in fractional

order Sobolev spaces for the p-Laplace equation on polyhedral domains. J.

Anal. Appl, 24:353-237, 2005.
[13] A. Elbert. A half-linear second order differential equation. In Qualitative

theory of differential equations, Vol. I, II (Szeged, 1979), volume 30 of Colloq.

Math. Soc. Jdnos Bolyai, pages 153-180. North-Holland, Amsterdam, 1981.

20

[14] L. C. Evans, M. Feldman, and R. F. Gariepy. Fast/slow diffusion and collapsing
sandpiles. J. Differential Equations, 137(l):166-209, 1997.

[15] M. Falcone and S. Finzi Vita. A finitc-diffcrcncc approximation of a two-
layer system for growing sandpiles. SIAM J. Sci. Comput., 28(3):1120-1132
(electronic), 2006.

[16] N. Igbida. A generalized collapsing sandpile model. Archiv der Mathematik,
94(2):193-200, 2010.

[17] A. Kuijpcr. Image analysis using p-Laplacian and geometrical PDEs. PAMM,

7(1):1011201-1011202, 2007.
[18] P. Lindqvist. Some remarkable sine and cosine functions. Ricerche Mat., 44

(2):269-290, 1995.

[19] W. Liu. On the stochastic p-Laplace equation. J. of Mathematical Analysis

and Applications, 360:737-751, 2009.
[20] M. Otani. A remark on certain nonlinear elliptic equations. Proc. Fac. Sci.

Tokai Univ., 19:23-28, 1984.
[21] C. Prevot and M. Rockner. A Concise Course on Stochastic Partial Differential

Equations. Springer, 2007.
[22] J. Rossi. Tug-of-war games. Games that PDE people like to play. Preprint,

2010.

[23] H-Y. Zhang, Q-C. Peng, and Y-D. Wu. Wavelet inpainting based onp-Laplace
operator. Acta Automatica Sinica, 33(5):546-549, 2007.

Department of Mathematics and Maxwell Institute for Mathematical Sciences
Heriot-Watt University, Edinburgh EH14 4AS, United Kingdom
E-mail address: 1 L . Boulton<3hw. ac .uk and 2 G . J . LordOhw. ac .uk